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A major obstacle to implementing Shor's quantum number-factoring algorithm is the large size 
of modular-exponentiation circuits. We reduce this bottleneck by customizing reversible circuits 
for modular multiplication to individual runs of Shor's algorithm. Our circuit-synthesis procedure 
exploits spectral properties of multiplication operators and constructs optimized circuits from the 
traces of the execution of an appropriate GCD algorithm. Empirically, gate counts are reduced by 
4-5 times, and circuit latency is reduced by larger factors. 
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I. INTRODUCTION 

Shor's number factoring remains the most striking al- 
gorithm for quantum computation as it quickly solves an 
important task [lj for which no conventional fast algo- 
rithms were found in 2,300 years [16J . Today, a scalable 
implementation of Shor's technique would have dire im- 
plications to Internet commerce. Laboratory demonstra- 
tions factored 15 = 3-5 circa 2000 [2], but further progress 
was slow [3-5 as factoring sizable semiprimes requires 
very large circuits. The bottleneck of Shor's number fac- 
toring is in modular exponentiation — a reversible cir- 
cuit computing (b z mod M) for known coprime integers 
b and M. This computation is performed as a sequence 
of conditional modular multiplications (CM) [3] by pre- 
compiled powers of a randomly selected base value (6), 
controlled by the bits of z (Figure [T]). In most cases, 6 = 2 
or b = 3 suffice [7 . Such CM blocks are assembled from 
unmodified unconditional modular multiplication blocks 
(1AM) using pre- and post-processing: since multiplica- 
tion always preserves the integer 0, a 1AM block can be 
"turned off" by conditionally swapping a with its in- 
puts and then restoring the inputs by an identical swap. 
Conditional swaps can be simplified, and further circuit 
optimizations focus on UM blocks. These steps are re- 
viewed in detail in [TJ. 

II. PRIOR WORK 

UM blocks are assembled from modular additions and 
multiplications by two, scheduled according to the bi- 
nary expansion of the constant multiplicand and its mod- 
ular inverse [B] (see a contemporary summary in [JJ). 
In one popular approach, the input value x is copied 
into a zero-initialized register to obtain (x,x). To com- 
pute 13ir, follow the binary expansion 13=61101: (x,x)- 
(2x,x)-(3x,x)-(6x,x)-(12x,x)-(13x,x). Now the second 
register must be restored to for the next UM block 
to use it. However, this requires dividing 13a; by 13, 
i.e., multiplying by the modular inverse of 13. For 
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M = 101113 = 569 • 1777, the inverse of C = 13 is 
77778, (IOOIOIIIIIIOIOOIO2), requiring a large circuit. 
In [7J, we constructed alternative circuits without com- 
puting modular inverses. To accomplish this, we intro- 
duced circuit blocks for modular multiplication and di- 
vision by two that restore their ancillae to 0. We then 
estimated costs of circuit blocks for modular addition, 
subtraction, multiplication and division by two, and sev- 
eral others [TJ Table 2]. Using these blocks, we found 
optimal UM circuits for each C, M up to 15 bits. The 
same procedure can be used for different cost estimates, 
but optimal search does not scale well beyond 15 bits. 
Numerical results demonstrated that traditional circuits 
based on binary expansion are far from optimal, thus 
asking for scalable constructions beyond 15 bits. 

Researchers optimizing circuits for Shor's algorithm 
[SI |9] adapted these circuits to use only nearest-neighbor 
quantum couplings [10] and restructured them to lever- 
age parallel processing . Applying multiple quan- 
tum couplings in parallel allows one to finish computa- 
tion faster. The smaller required lifespan of individual 
qubits additionally reduces the susceptibility of qubits 
to decoherence and decreases the overall need for quan- 
tum error-correction. The runtime (latency) of parallel 
quantum computation is estimated by the depth of its 
quantum circuit, i.e., the maximum number of gates on 
any input-output path. Depth reductions in the litera- 
ture sharply increase the required number of qubits, e.g., 
by 50 times or more, making them impractical for mod- 
ern experimental environments where controlling 50-100 
qubits remains a challenge. Vice versa, prior circuits with 
1-2 fewer qubits use more gates |12j . Rosenbaum has 
shown |13| how to adapt unrestricted circuits to nearest- 
neighbor architecture using teleportation, while asymp- 
totically preserving their depth. For Shor's algorithm, 
the range of practical interest is currently between several 
hundred and a thousand logical qubits, where FFT-based 
multiplication needs more gates than simpler techniques. 

Our circuits moderately increase qubit counts to sig- 
nificantly decrease gate counts and circuit depth. Built 
from standard components, they are readily adapted to 
nearest-neighbor quantum architectures by optimizing 
these components to each particular architecture [10J. 
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Figure 1: (a) Modular exponentiation using conditional modular multiplications [6,. (b) Conditional multiplications imple- 
mented using unmodified unconditional modular multiplication blocks and conditional swaps with a zero register [7J. 



III. NEW CIRCUITS 

We propose two-register UM circuits to compute (Cx 
mod M) for coprime C and M, GCD(C,M) = 1. These 
circuits transform (a;, 0) into (x,x) using parallel CNOT 
gates and then compute (Cx mod M, 0). Clearing ancil- 
lae in the second register allows the next circuit module 
to use them again. As reversible building blocks, we use 
modular addition and subtraction between the two reg- 
isters (a, b) > (a ± b mod M, b) and (a, b) i— > (a, a ± b 
mod M), as well as circuits from [7] for modular multipli- 
cation by two that clear their ancillae a H> 2a mod M. 

Our key insight is to use the coprimality of C and 
M (guaranteed in Shor's algorithm) to read off a cir- 
cuit from the execution trace of an appropriate GCD 
algorithm. Recall that the Euclidean GCD algorithm 
for (A, B) proceeds by replacing the larger number A 
with (^4 mod B) until the result evaluates to 0. For 
C = 13 and M = 21, this produces the Fibonacci se- 
quence (21,13)-(8,13)-(8,5)-(3,5)-(3,2)-(l,2)-(l,0). For 
convenience, one may consider the last configuration to 
be (1,1), so that each step performs a subtraction — a 
simpler operation than modular reduction. As a result, 
we obtain GCD(C, M) = 1. Reversing the order of opera- 
tions, interpreting each number as a multiple of x starting 
with (lx, lx), and mapping each step into a mod-21 ad- 
dition, we obtain (x, x)-(2x, x)-(2x, 3x)-(5x, 3cc)-(5a;, 8a;)- 
(13a:, 8a;)-(13x, 21a;) = (13a;, 0). Since 21x mod 21 = 0, 
the second register is restored to 0. This UM circuit by- 
passes Bennett's construction based on modular inverses 
(Section [II]) and is smaller than prior art [1, 6J. 

Unfortunately, some modular reductions in the Eu- 
clidean GCD algorithm may require a large number of 
gates. Consider (11a; mod 21) and its Euclidean GCD 
trace (21, 11)-(10, 11)-(10, 1)-(1, 1). Implementing the 
last mod operation by nine subtractions produces a 
sequence of nine mod-21 additions (x, x)-(2x, x)-(3x, x)- 
. . .-(10x, x). To improve efficiency, we replace the Eu- 
clidean GCD algorithm by a binary GCD algorithm that 
avoids the mod operation and uses a shortcut for the 
case of odd GCD. Given a pair of odd numbers, the larger 
one is replaced by their difference, which must be even. 
Any even number is divided by two, which can be imple- 
mented by a controlled bit-shift (as shown in [7]). 
For even A and B, (A, B) = (A/2, B/2) 
For even A and odd B, (A, B) = (A/2, B) 
For odd A and even B, (A, B) = (A, B/2) 
For odd A and B, if A < B, (A, B) = (A,B- A) 
else (A, B) — (A — B, B) 



One stops when A=B=GCD=1 (assuming coprime in- 
puts). The sequence of operations performed for our 
example (21, 11)-(10, ll)-(5, ll)-(5, 6)-(5, 3)-(2, 3)-(l, 3)- 
(1, 2)-(l, 1) can be improved by (2, 3)-(2, 1)-(1, 1). To ob- 
tain a circuit, such sequences are reversed and interpreted 
as modular multiplications by two and modular addi- 
tions, with the initial state (la;, la;). Further improve- 
ments are obtained by allowing both subtractions and ad- 
ditions, e.g., 15a; = 16a; — x versus 15a; = 8x+4a; + 2a; + l3: 
(here 16a;, 8a;, etc are computed by doubling). 

In a more involved example (7x mod 1017), the ad- 
dition leading to (7, 1024) is a better first step than the 
subtraction leading to (7, 1010), because (7, 1024) enables 
eight successive divisions by two which reduce the values 
down to (7, 4) faster than subtractions would. Then, sub- 
tractions become the best operators: (3, 4)-(3, l)-(2, 1)- 
(1,1). This optimization relies on a three-step lookahead. 
To select each next operator, we consider all possible irre- 
dundant three-step sequences of operators (modular ad- 
dition, subtraction and division by two), find their final 
states, and score the remaining circuit according to the 
trace of the binary GCD algorithm (without lookahead) . 
The cost of each operator/step can be specific to the 
quantum machine. Taking the best three-step sequence, 
we commit to its first operator. The remaining two steps 
are ignored, and the next operator is chosen by a sepa- 
rate round of lookahead. For (11a; mod 21), we obtain 
(21, 11)-(10, ll)-(5, ll)-(5, 6)-(5, l)-(4, l)-(2, 1)-(1, 1). 

IV. EMPIRICAL VALIDATION 

Our algorithms for on-demand construction of modu- 
lar multiplication circuits |17j were embedded into the 
framework of Figure ^ The number of ancillae in re- 
sulting mod-exp circuits was 5n + 2 (as in (7|), but sev- 
eral optimizations from [7] were not used, and the num- 
ber of mod- mult blocks was exactly as in j3] . Our soft- 
ware was written in C++ using the GNU MP library 
(for multi-precision arithmetic) supplied with the GCC 
4.6.3 compiler on Linux. We used a workstation with an 
Intel® Core™ 2 Duo 2.2 GHz CPU and 2 GB Memory. 

To evaluate our optimizations of Shor's number- 
factoring algorithm, we studied all odd n-bit semiprimc 
values of the modulus (M = pq) for 7 < n < 15, and a 
subset of n-bit M values for n =16-512 that are prod- 
ucts of the 1st and 10th largest n/2-bit primes. Circuit 
sizes for n < 16 were averaged over all M-coprime C 
values. Results for n = 16 include all coprime C values 
for the given M. For 24-, 32-, 48-, and 64-bit M val- 



3 



Table I: Circuits produced by our technique and prior art, compared by Toffoli gate counts. Circuit sizes for n < 16 are averaged 
over all Af-coprime C values. Results for n = 16 include all coprime C values for the given M. For 24-, 32-, 48-, and 64-bit 
M values, results are averaged over the first 5000 coprime C values. For larger n values in mod-mult circuits, only C — —1/17 
mod M (e.g., 47679095568306588235294117647058823529411764705882352941176470588235294117647 for n = 256) are shown. 
For modular exponentiation, results include all C values appearing in UM blocks for b = 2. All results reported are circuit 
sizes (Toffoli gate counts), except for values in the depth column. For 'Avg ratio' in mod-mult, we used Ours/[7] and [6]/Ours. 
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ues, results were averaged over the first 5000 coprime C 
values. For larger n values in modular multiplication cir- 
cuits, only C = —1/17 mod M are shown. Results for 
modular exponentiation include all C values appearing 
in unconditional modular multiplication blocks for 6 = 2 
(Figure These are C = b 2 °%M, C = b 2± %M, 
and C = 6 22 " _1 %M. For n < 15, Table [j] shows that 
circuits found by our heuristic are closer to optimal cir- 
cuits [7J than to scalable circuits from [5j. Beyond the 
reach of optimal techniques (n > 24) , Figure [2] shows 
that our circuits are at least 4.5 times smaller and retain 
their advantage as n increases. Our runtimes ranged from 
negligible (n < 32) to 30 min for one 512-bit (M, C) pair. 

To compare our circuits with latency(depth)-optimized 
constructions in |11| , we note that the most accurate data 
in fllj are given for n = 128. Our smallest 128-bit mod- 
exp circuits use 1.97 x 10 7 Toffoli gates with 642 ancillae. 
To reduce the latency of our circuits, we replaced linear- 
depth Cuccaro adders with log n-depth adders from [2] 
also used in [TT]. Accordingly, circuit depth is reduced to 
2.03 x 10 6 Toffoli gates with - 900 ancillae. This process 
is outlined in the next section, but here we summarize 
the results. A circuit with 660 ancillae |lll Algorithm G, 
Table II] exhibits latency 1.50 x 10 7 Toffoli gates. The 
best circuit in [11, Algorithm E, Table II] has latency 
1.71 x 10 5 Toffoli gates but uses 12657 ancillae, which 
is far less practical with technology under development 
today. Circuit depths of our modular exponentiation cir- 
cuits for all attempted n values are reported in Table |T] A 
quantum machine with only some limited form of paral- 
lelism may still benefit from our techniques, given strong 
results for both parallel and sequential cases. 



V. REDUCING CIRCUIT LATENCY 

Our circuits can be adapted to quantum architec- 
tures with high degree of parallelism by replacing build- 
ing blocks by parallelized variants. Circuit-size calcula- 
tions in Table |T] are based on the costs of circuit mod- 
ules (addition, subtraction, modular multiplication by 2, 
etc) from [71 Table 2]. Cuccaro adders used in [7] are 
small, but exhibit linear latency. To optimize latency 
for comparisons to [TT], we replaced Cuccaro adders with 
QCLA adders from [2] (also used in [TT]) whose depth 
is (41og 2 n + 3, 4, 2) in terms of (T,C,N) gates. As in 
[TT], we measure latency in Toffoli gates. This results 
in circuit latency (depth) 4 log 2 n + 3 for additive oper- 
ators (~l,~2,+l,+2,-l,-2) i n [7J Table 2]. The oper- 
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Figure 2: Asymptotic behavior of circuit-size ratios between 
Beckman et al |H| and our constructions. 
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ators that perform modular multiplication (dl,d2) and 
division by two (hi , h2) exhibit latency 6 log 2 n + 12. To 
count the number of ancillae in our modular exponen- 
tiation circuits, note that QCLA adders from [T3] need 
2n — \ogn — 2 ancillae (vs. 1 for Cuccaro adders). Given 
that QCLA adders clear all ancillae, the number of an- 
cillae in our mod-exp circuits grows to < In. 

We also restructure n one-bit controlled-SWAP gates 
with shared control to reduce latency from n to log 2 n. 
The control bit is temporarily copied to n zero-initialized 
ancillae (with log 2 n latency) [IS] . We use n parallel one- 
bit controlled-SWAP gates, and then clear the n ancilla 
(also with log 2 n latency) . Because these ancillae are 
cleared immediately, we can share them with the QCLA 
ancillae. Thus, the overhead is 21og 2 n latency and 2n 
CNOT gates used to copy/clear ancillae. 

VI. CONCLUSIONS 

The n-bit multiplication circuits developed in this work 
significantly simplify the implementation of Shor's algo- 
rithm, but use 8(n 2 ) gates, as do traditional circuits. 
Circuit sizes are improved by large constant factors. 
These factors appear exaggerated for small qubit arrays 
because, for 6 = 2, our construction implements a non- 
trivial fraction of modular multiplications using 0(n) 
gates using circuit blocks from j7j. In contrast, prior 
work typically uses generic 0(n 2 ) circuits regardless of b. 
We have experimented with several enhancements to our 
technique, but the resulting improvement was not justi- 



fied by the increased runtime and programming difficulty. 

Connections between number factoring and CCD com- 
putation were known to Euclid around 300 B.C.E. Today 
the two problems play similar roles in their respective 
complexity classes. Number- factoring is in NP (prob- 
lems whose solutions can be checked in polynomial time), 
not known to be in P (problems solvable in polynomial 
time), but is not believed to be NP-complete (most dif- 
ficult problems in NP). GCD is in P, not known to 
be in NC (problems that can be solved very efficiently 
when many parallel processors are available), but is not 
believed to be P-complete (inherently sequential). Un- 
like provably-hard problems (such as Boolean satisfiabil- 
ity), or problems for which fast serial and parallel algo- 
rithms are known (such as sorting) , number factoring and 
GCD appear to be good candidates for demonstrations 
of physics-based computing that exploits parallelism. 

Our use of GCD algorithms to speed up modular expo- 
nentiation and number factoring incurs only small over- 
head. All the invocations of our GCD-based circuit con- 
struction for one run of Shor's algorithm can run in par- 
allel because they are independent. Thus, the classical- 
computing overhead of our technique for one run of Shor's 
algorithm is limited to 1-2 GCD-based circuit construc- 
tions. This overhead is acceptable because Shor's al- 
gorithm performs multiple CCD-like computations after 
quantum measurement. 
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